rm(list=ls())
library(foreign)
library(doMC)
library(mvtnorm)
require(foreach)
require(doMC)
registerDoMC(cores=3)
library(rjags); load.module("glm"); load.module("lecuyer")

set.seed(100)

statadata <- read.dta("2_forjags.dta")
summary(statadata)
nrow(statadata)

# remove 1952 (1st quarter 1952 missing)
statadata <- statadata[4:nrow(statadata),]

## Re-scale for use on 0-1 scale of macro
statadata$partyics <- statadata$partyics/100
statadata$partyapprove <- statadata$partyapprove/100
statadata$lagpartyics <- statadata$lagpartyics/100
statadata$lagpartyapprove <- statadata$lagpartyapprove/100

ols.fit <- lm(adjusted_gallup~lagmacro + party + partyics + partyapprove + lagparty + lagpartyics + lagpartyapprove,data=statadata)
summary(ols.fit)


b0 <- coef(ols.fit)[-1]
B0 <- diag(1E-04,7)
phi.sig <- diag(7)*.1

forJags <- list(
  y = statadata$adjusted_gallup,
  ylag = statadata$lagmacro,
  party = statadata$party,
  partyics = statadata$partyics,
  partyapprove = statadata$partyapprove,
  lagparty = statadata$lagparty,
  lagpartyics = statadata$lagpartyics,
  lagpartyapprove = statadata$lagpartyapprove,
  macrovar =statadata$macrovar,
  T = nrow(statadata),
  b0 = b0,
  B0 = B0
  )


seeds<- c(1001,1002,1003)
jinits <- function(i){
    list(.RNG.name="lecuyer::RngStream",
         .RNG.seed=seeds[i],
         "munaught" = runif(1,.4,.6),
         "phi" = as.vector(rmvnorm(1,mean=b0,sigma=phi.sig)),
         "sigma.a" = runif(1,.008,.016),
         "sigma.b" = runif(1,.01,.02),
         "epsilon" = runif(1,0,.5)
         )
}
 

jags.parsamples <- foreach(i=1:getDoParWorkers(),.combine="c") %dopar% {
   
  foo <- jags.model(file="mccullochtsay_fulladl.jags",
                  data=forJags,
                  n.adapt=20000,
                  n.chains=1,
                  inits=jinits(i))
  update(foo,n.iter=80000,by=100,progress.bar="text")
  out <- coda.samples(foo,
                    variable.names=c("sigma.a","sigma.b","epsilon","phi","mu","delta","yexp"),
                    n.iter=2000000,
                    thin=500,
                    by=1000,
                    progress.bar="text")
  return(codaResults=out)
}

class(jags.parsamples) <- "mcmc.list"

summary(jags.parsamples)

library(coda)
pdf("TableC3_MSAR_ADL.pdf")
plot(jags.parsamples[,'sigma.a'],main="Sigma Alpha")
plot(jags.parsamples[,'sigma.b'],main="Sigma Beta")
plot(jags.parsamples[,'phi[1]'],main="Lag Macro")
plot(jags.parsamples[,'phi[3]'],main="ICS")
plot(jags.parsamples[,'phi[4]'],main="Pres Approval")
plot(jags.parsamples[,'phi[6]'],main="ICS Lag")
plot(jags.parsamples[,'phi[7]'],main="Pres Approval Lag")
plot(jags.parsamples[,'epsilon'],main="Epsilon")
plot(jags.parsamples[,'mu[120]'],main="Mu 120")
dev.off()

save.image(file="TableC3_MSAR_ADL.Rdata")

rm(list=ls())

                   

